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Summary 

A viscoplastic stress-strain analysis of an experimental 
cylindrical thrust chamber is presented. A viscoplastic consti- 
tutive model incorporating a single internal state variable that 
represents kinematic hardening was employed to investigate 
whether such a viscoplastic model could predict the experi- 
mentally observed behavior of the thrust chamber. Two types 
of loading cycles were considered: a short cycle of 3.5-sec 
duration that corresponded to the experiments, and an extended 
loading cycle of 485. 1 -sec duration that is typical of the Space 
Shuttle Main Engine (SSME) operating cycle. The analysis 
qualitatively replicated the deformation behavior of the 
component as observed in experiments at NASA Lewis 
Research Center. These experiments were designed to simulate 
the SSME operating conditions. The analysis also showed that 
the mode and location of failure in the component may depend 
on the loading cycle. The results indicate that using viscoplastic 
models for structural analysis can lead to a more realistic life 
assessment of thrust chambers. 

Introduction 

Progressive deformation and thinning of coolant-channel 
walls in high-pressure, reusable rocket thrust chambers have 
been observed during operation of engines such as the Space 
Shuttle Main Engine (SSME). This inelastic ratcheting behavior 
is attributed to the significant, thermally induced inelastic 
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strains that occur in the hot-gas-side wall during cyclic firing, 
as well as to the biasing pressure differential across the wall. 
After numerous thermal cycles, ratcheting strains accumulate 
to the extent that cracks form (i.e., failure occurs) in the 
cooling passage wall. Numerous analytical and experimental 
studies have been undertaken in recent years to examine this 
phenomenon (refs. 1 to 5). 

For an accurate estimation of the life of thrust chambers, 
a realistic inelastic stress-strain analysis is required. Conven- 
tional inelastic analyses of components treat the plastic (time- 
independent) and creep (time-dependent) strains as independent 
noninteracting entities; these are inadequate for this appli- 
cation, since at elevated temperatures, interactions are known 
to exist between plastic and creep strains within metallic 
materials (ref. 6). Alternatively, unified viscoplastic analyses 
provide realistic descriptions of high-temperature inelastic 
behavior of materials, in that all inelastic strains (e.g., creep, 
plastic, relaxation, and their interactions) are accounted for 
as a single, time-dependent quantity. 

In the present study a unified viscoplastic finite element 
stress-strain analysis was made of an experimental thrust 
chamber composed of NARloy-Z, a copper-based alloy. The 
objective of this study was twofold: (1) to qualitatively predict 
the experimentally verified thinning of the cooling passage and 
the so-called “dog house” effect (see ref. 2) by using a 
viscoplastic model put forth by Robinson (ref. 7); and (2) to 
investigate the influence of loading cycle duration on the 
structural response of the cylindrical chamber. To achieve the 
second objective, two types of thermomechanical loading 
cycles (see table I) were employed— one corresponding to the 
experimental simulations and one being more typical of the 
cycle seen by the SSME. 
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TABLE I. -CYCLIC TEMPERATURE AND PRESSURE 
LOADING HISTORIES 

[From Ref. 3.] 


(a) Duration of Regions 


Region 

Cycle time, sec 


Short 

Extended 

I 

0.10 

22.50 

II 

.90 

.90 

III 

.80 

180.00 

IV 

.45 

.45 

V 

1.25 

281.25 

Total cycle 



time 

3.50 

485.10 


a 

Heat transfer 
coefficient, 

temperature, phase 
or pressure I 


Typical Cycle History 
m 



Time, sec 


(b) Cycle Phase 



Cold 

Hot 

Hot-gas-side heat transfer coefficient. 



W/cm 2 -K (Btu/in. 2 -sec-R) 

0 

2.02 (0.00685) 

Hot-gas-side adiabatic wall 
temperature, K (°R) 

278 (500) 

3364 (6055) 

Hot-gas-side wall pressure, 
kN/m 2 (psia) 

Coolant-side heat transfer coefficient, 

96.5 (14.0) 

2780 (403) 

W/cm 2 -K (Btu/in. 2 -sec-R) 

10.2 (0.0345) 

4.83 (0.0164) 

Coolant-side bulk temperature, K (°R) 

28 (50) 

50 (90) 

Coolant-side wall pressure, 
kN/m 2 (psia) 

1 5100 (740) 

6550 (950) 


A cross section of the experimental cylinder wall, which 
contains 72 cooling channels, is shown in figure 1 . The shaded 
area in the figure depicts the segment of the component that 
was modeled to perform the thermal and structural analyses. 
The cyclic thermal and pressure loadings (ref. 3) shown in 
table I lead progressively to a deformed segment that results 
in the “dog house” effect noted in reference 2. Figure 2, which 
shows this dog house, illustrates how the component fails by 
thinning of the cooling channel wall and eventual tensile 
rupture. Further details of experiments and apparatus can be 
found in references 2 and 3. 

This paper begins with a brief description of Robinson’s 
viscoplastic model and the finite element analysis that was 
performed. Results are then presented and discussed, and 
conclusions are drawn. Future work is suggested. 



Figure 1. — Cylinder wall cross section showing analyzed segment and 
dimensions. 



Figure 2.— Cross section of 1 /2-hard Amzirc cylinder at the throat plane after 
393 thermal cycles. 


Robinson’s Viscoplastic Model 

Robinson’s model (ref. 7) employs a dissipation potential 
to derive the flow and evolutionary laws for the inelastic strain 
and internal state variables. The model incorporates a single 
internal state variable representing kinematic hardening. The 
material behavior is elastic for all the stress states within the 
dissipation potential and is viscoplastic for the stress states 
outside. A small displacement and a small strain formulation 
are employed. 

The total strain rate is decomposed into elastic ifj, 
inelastic (including plastic, creep, relaxation, etc.), and 
thermal ef strain-rate components. Thus 

e,j = (i,j = 1 , 2 , 3 ) ( 1 ) 

The elastic strain rate for an isotropic material is governed 
by Hooke’s law: 


-el _ 1 + v . v . . 

e ij ~ °ij ~ £ °kk O/j 


( 2 ) 
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where E is the Young’s modulus, v is the Poisson’s ratio, and 
a i} is the stress. The repeated subscripts in equation (2) and 
elsewhere imply summation over their range, and h i} is the 
Kronecker delta function. A dot over a symbol denotes its 
derivative with respect to time t. 

The nonisothermal multiaxial inelastic constitutive equations 
for the model are given below: 


Flow Law 


AF* Xij 
V7? 


F> 0 and Sy E/, > 0 


0 F < 0 or F < 0 and S {J E f> * < 0 


In the preceding equation, is the deviatoric stress, and a {j 
and K are the internal state variables. The variable a ijy called 
the back stress, accounts for the kinematic hardening, whereas 
K , taken here to be a constant, is called the drag stress and 
represents the isotropic hardening of the material. The value 
of K at the reference temperature is denoted by K n , and the 
minimum value attainable by G is denoted by G a . The 
inequalities in equations (3) and (4) define boundaries across 
which the flow and evolutionary laws change form discon- 
tinuously. To facilitate the numerical computations, these 
boundaries are smoothed by using a spline function given in 
the next section. Seven numerical parameters are required to 
characterize a given material. These parameters are A t n f m, 
H f R, G a , and K(K a ). The values of these parameters for the 
copper alloy NARloy-Z, taken from reference 8, are listed 
in table II. 


(3) 


Evolutionary Law 

m H RG m ~ 0 


^ w V7? y 


H RG'r 0 

; e ? a, 


G < G a and 5,y a XJ > 0 


G < G n 


Lpo U vT 2 lJ or G > G a and S tj a < 0 


where 


(4) 


f = - 2 - 1 

K 2 


G = 


Ki 


J 2 ~ ~ E y 


h — - a ij a ij 


(5) 

(6) 

(7) 

(8) 


TABLE II — VALUES OF MATERIAL 
PARAMETERS" FOR ROBINSON'S MODEL 


[From ref. 8.] 


Parameter 

Value 

A 

1 .385 x 10 -8 

n 

4 

m 

4.365 

0 

0.533 x 1CT 6 T 2 + 0.8 

Qo 

40 000 

T n 

811 K (1000 °F) 

G„ 

0.04 

fO- 

69.88 - 0.067 T 

AV 

K(T„) 

H 

(1.67x 10 4 )(6.895) ,5+ '/(3 Kj) 

R 

(2.19X 10 _, )(6.895) ,+ ^ _m (3 


xexp |q o 

E 

1.47x10 s - 70.5 XT 

V 

0.34 


a Unils are as follows; T is in Kelvin; n, m, 0, and G () arc unitless; and 
the other parameters are consistent with the units of MPa and hr 


— Sy a :j 


(9) 


Sjj Ojj 5jj 

a ij ~ a ij ~ ~ °kk ?>ij 
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Spline Function 

The discontinuous boundaries in Robinson’s viscoplastic 
(10) model should be smoothed to facilitate numerical computa- 
tions. This smoothing is achieved by defining a “spline 
function” P{x) (refs. 7 and 9) on the interval (-1,1) as 
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-1 < .* < 0 


Backside wall 


r 

(1 + *) 2 


P(X) = < 


1 - 


1 


(1 -X) 2 


0 

V. 


0 < x < 1 


( 11 ) 


X > 1 
* < -1 


Analyses 

The details of the finite element model, and thermal and 
stress analyses are provided in this section. 

Finite Element Model 

Figure 3 depicts the finite element model used for the thermal 
and stress analyses. It consists of 35 elements and 54 nodes. 
To facilitate a quasi-three-dimensional analysis, generalized 
plane strain elements were used to model the smallest repeating 
segment of the cylinder wall. Because the wall was symmetrical, 
only one-half of a cooling channel was modeled. 

Thermal Analysis 

In reference 2, the thermal analyzer program SIND A (ref. 10) 
was used to perform the thermal analysis and to generate the 
time-dependent cross-sectional temperatures. Time-dependent 
hot-gas-side and coolant-side boundary conditions were input 
to the SINDA program to obtain the temperature profile as 
a function of time. To estimate the quantities for which no 
experimental data were available, some assumptions and 
approximations were made in the thermal analysis. Complete 
details of the thermal analysis are provided in reference 2. 

Finite Element Analysis 

By employing the finite element program MARC (ref. 1 1) 
and a numerical solution procedure put forth by Arya et al. 
(refs. 12 and 13), a stress-strain analysis of the cylindrical 
thrust chamber was performed. A self-adaptive integration 
strategy, based on the explicit forward Euler method and 
described in reference 14, was employed for the time integra- 
tion of the nonlinear and mathematically “stiff’ constitutive 
equations of the viscoplastic model. The integration strategy 
is easy to implement, and it has been employed successfully 
for the solution of problems involving complex geometries and 
loadings (refs. 4 and 15). 

The finite element thermostructural analysis was performed 
with the time-dependent temperature distribution obtained from 
the thermal analysis described in the preceding section. Time- 


Cooling 

channel 



Figure 3.— Finite element model for cylindrical thrust chamber; 35 elements, 
54 nodes. 


varying nodal temperatures and elemental pressure loadings 
were applied to trace the loading cycles of table I. Two loading 
cycles of different durations were used for numerical compu- 
tations. The short cycle of 3.5-sec duration corresponded to 
the loading cycles of the experiments. The extended loading 
cycle of 485.1 -sec duration was typical of SSME operation. 
The latter cycle offered the opportunity to study the effects 
of time-dependent creep deformation on the structural response 
of the cylindrical chamber liner. The nonlinear variations in 
the temperature-dependent material properties were accounted 
for in the computations. The temperature distribution obtained 
from the thermal analysis described in the preceding section 
was used for the stress analysis. 
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Results and Discussion 

Short Loading Cycle 

The circumferential (x-direction) stress distribution in the 
segment of the component used for the finite element analysis 
is shown in figure 4 at the completion of the 1st and 12th short 
loading cycles. The color figures included in this paper are 
aimed at providing a quick inference with regard to the tensile 
or compressive nature of the stress (or strain) distribution in 
the segment. The stress (or strain) distributions shown in the 
color figures thus represent the segment's qualitative response 
to cyclic thermomechanical loading rather than its quantitative 
response. The circumferential stress is tensile and, in general, 
decreases with the number of loading cycles. The radial 
(y-direction) stress distribution is shown in figure 5; here some 
redistribution of the radial stress with loading cycles is observed 
in the upper part of the segment. In the lower part of the 


segment the radial stress remains unaltered at some points, 
but it increases in magnitude at other points of the segment; 
for example, in the vicinity of point B. 

The circumferential (x-direction) and radial (y-direction) 
mechanical strain distributions in the segment at the end of 
the 1st and 12th short loading cycles are shown in figures 6 
and 7, respectively. The circumferential strains are compressive 
everywhere and are found to increase in magnitude in the 
neighborhood of points A and D with the number of loading 
cycles. The radial (y-direction) strains are compressive in the 
segment at the end of the first cycle. Comparison of the radial 
strain distributions in the segment at the end of the 1st and 
12th loading cycles (see fig. 7) shows that the radial strain 
becomes less compressive with the loading cycles. 

The predicted shapes of the segment after the 1st and 12th 
short loading cycles are shown in figure 8. These figures are 
magnified by a factor of 50 so that the small difference in 
deformity can be seen. 


Stress, 
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300 


ksi 

43.5 



Figure 4.— Short cycle circumferential (.v-direction) stress distribution in segment. 
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Figure 5.— Short cycle radial (y-direction) stress distribution in segment. 
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Figure 6 — Short cycle circumferential (^-direction) strain distribution in segment. 
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Figure 7.— Short cycle radial (y-direction) strain distribution in segment. 
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Figure 9.— Extended cycle circumferential (.r-direction) stress distribution in segment. 


Extended Loading Cycle 


The circumferential (x-direction) and radial (v-direction) 
stress distributions in the segment at the end of the 1st and 
12th extended loading cycles are shown in figures 9 and 10, 
respectively. The circumferential stress is tensile and decreases 
with the loading cycles. The radial stress either remains 
unaltered or increases in magnitude with the loading cycles. 
A comparison of the radial and circumferential stress distri- 
butions for the short (figs. 4 and 5) and the extended (figs. 9 
and 10) loading cycles shows that, in general, the stress values 
for the extended loading cycles are lower in magnitude than 
the corresponding values for the short loading cycles. 

Figures 1 1 and 12 portray the circumferential (x-direction) 
and radial (y-direction) strain distributions in the segment at 
the completion of the 1st and 12th extended loading cycles. 


The circumferential strain increases in magnitude in the lower 
right part of the segment in the neighborhood of point D, 
whereas it decreases in magnitude in the lower left part of the 
segment in the vicinity of point B and increases in the neighbor- 
hood of point A. The radial strain becomes less compressive 
with the loading cycles at all points of the segment. The radial 
strain in the neighborhood of point D changes to tensile 
(12th cycle) from compressive (1st cycle). A comparison of 
the strain distribution for the 12th short (fig. 6) and extended 
(fig. 11) loading cycles reveals that, for the extended cycle, 
the magnitude of the circumferential strain is lower in the 
vicinity of points A and B and higher in the vicinity of points 
C and D than it is for the short cycle. The radial strain for 
the 12th extended loading cycle (fig. 12) is found to be the 
same or lower in magnitude than that for the short loading 
cycle (fig. 7) at the corresponding points of the segment— 
except near point B, where the trend is reversed. 
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Figure 10.— Extended cycle radial ^-direction) stress distribution in segment. 
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Figure 11. — Extended cycle circumferential (.v-direction) strain distribution in segment. 
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Figure 12.— Extended cycle radial ( y -direction) strain distribution in segment. 
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1st Cycle 12th Cycle 

Figure 13.— Predicted shapes of segment after different cycles (extended cycle; magnification factor * 50). 


Figure 13 depicts the deformed shapes of the segment after 
the 1st and 12th extended loading cycles. This figure shows 
noticeable deformation of the coolant-channel wall after the 
12th extended loading cycle. 

Figure 14 displays the circumferential (jc-direction) dis- 
placements at point C of the segment (see inset) as a function 
of the number of cycles for both the short and extended loading 
cycles. The displacement values plotted are the values at the 
completion of the corresponding loading cycle. Since point 
D of the segment (which lies in the same horizontal plane as 
the point C) is constrained in this direction, the values shown 
in figure 14 also give the change in thickness between points 
C and D. A comparison of the curves for the short and extended 
cycles reveals that, for the extended cycle, there is a greater 
thinning across the section CD of the segment. This fact is 
supported by examination of the circumferential strain ranges 
across the section CD for the short and extended loading 
cycles. Since the strain range for the extended cycle is larger 
(and compressive), greater thinning of the channel wall across 
section CD is implied for the extended cycle. 

The radial (y-direction) displacements at two locations of 
the segment (designated by points A and B in the inset) are 
shown in figure 15. The displacement values are plotted with 
respect to the number of cycles. Note that the displacement 


values at both locations A and B are higher in magnitude for 
the extended loading cycle. This means that the coolant-channel 
wall bulges out more for the extended cycle than for the short 
loading cycle. The difference in the displacement values at 
points A and B is the change in wall thickness (thinning) across 
the section AB of the segment. The absolute value of this 
difference, denoted by Ad , is plotted in figure 16 as a function 


1.5X10 -4 



Figure 14.— Decrease in circumferential thickness of segment wall at point C 
for both short and extended cycles. 
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Figure 15. — Radial displacement of points A and B. 
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Figure 16. — Decrease in radial thickness across plane AB of segment, for 

both short and extended cycles. 

of the number of cycles. It is interesting to note that the change 
in wall thickness of the coolant channel (thinning) for the 
extended loading cycle is smaller than that for the short loading 
cycle. Close examination of figures 7 and 12 reveals that at 
the end of 12th cycle the radial strain range across the section 
AB of the segment is smaller for the extended cycle than it 
is for the short cycle. The strains are compressive. This smaller 
(compressive) radial strain range indicates a smaller thinning 
of the coolant-channel wall per cycle for the extended cycle 
history. Since, for the extended cycle, (1) the displacement 
values are higher in the radial direction and (2) the wall 
thinning across section AB is smaller in the radial direction 
and larger across section CD in the circumferential direction 
than that for the short loading cycle, the failure modes and 
their location in the segment may be different for the two types 
of loading cycles. 

The circumferential fv-direction) stress distribution in the 
segment at the hot phase of the 1st and 12th extended loading 
cycles is shown in figure 17. Figure 18 depicts the corre- 
sponding temperature distribution. We can see from figure 17 


that, at the hot phase of the 12th extended cycle, the maximum 
(in magnitude) circumferential stress is compressive and occurs 
across section AB. The circumferential stress at the completion 
of the same cycle (fig. 9) is also at a maximum across section 
AB, but it is tensile. Since a compressive stress cannot cause 
the radial cracks seen in the experiments (fig. 3), the cracks 
(and hence the failure of the segment) are unlikely to occur 
at the hot end of the cycle. The cracks (and failure) may occur 
only at the end of the cycle or at the beginning of the next 
cycle, since the circumferential stress then becomes tensile. 

Figure 19 shows the predicted shapes of the NARloy-Z 
segment at the end of 24th short and the 4th extended loading 
cycles. By employing Robinson’s viscoplastic model, we were 
able to predict the “dog house" effect observed in the experi- 
ments of reference 2 (see fig. 2). It is important to note that 
Robinson’s model has only one internal state variable. In an 
earlier analysis of a copper cylindrical thrust chamber, Arya 
(ref. 4) employed a viscoplastic model developed by Freed 
(ref. 16) that incorporated two internal state variables, rep- 
resenting kinematic and isotropic hardening, to predict the 
“dog house" effect. These results indicate that for the present 
problem a viscoplastic model with a single internal state 
variable, representing kinematic hardening, can be successfully 
employed to predict the experimentally observed behavior. In 
fact, for a nonisothermal loading cycle such as that used in 
this work, the isotropic hardening variable in Freed’s model 
saturates quickly; this explains why the two models predict 
similar deformation behavior of the segment. The predicted 
shape for the copper thrust chamber at the end of the fourth 
short loading cycle is also shown in figure 19. This figure 
shows a larger deformation for the copper chamber than for 
the NARloy-Z chamber, as would be expected since copper 
is a softer material than NARloy-Z. 

Thermal Loading Cycle 

To assess the role that pressure loading plays on the “dog 
house" effect, the segment was analyzed by considering 
thermal loading only. The predicted shapes of the segment at 
the end of the 1st and 14th short loading cycles are shown 
in figure 20. In this figure the segment bulges inward when 
there is no pressure loading. This implies that the “dog house" 
effect seen in the experiments and qualitatively predicted from 
the current viscoplastic analysis is due to the presence ot a 
biasing pressure differential across the coolant-channel wall. 
Although not shown here, intuitively similar conclusions may 
be drawn for the extended loading cycle also. 

Research Work in Process 

Owing to their light weight and enhanced strength, metal- 
matrix composite (MMC) materials are becoming increasingly 
important in aerospace applications. Currently, an MMC 
cylindrical thrust chamber is being analyzed by using the 
transversely isotropic viscoplastic model put forth by Robinson 
(D.N. Robinson, 1989, University of Akron, Akron, Ohio, 


15 




MPa 

50 


Stress, 


ksi 

7.3 




Figure 17.— Extended cycle circumferential fx-di recti on) stress distribution in segment, for hot phase of cycle. 
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Figure 20.— Predicted shapes of segments owing to thermal loading only (short cycle). 


personal communication). Elastic contribution in the model 
is being calculated by using a transversely isotropic thermo- 
elastic theory put forth by Arnold (ref. 17). The study reported 
herein was partly motivated by the current MMC thrust 
chamber analysis since the transversely isotropic model, which 
is similar to Robinson's isotropic model, incorporates only a 
single internal state variable. Thus, before analyzing the MMC 
chamber, we first had to ascertain that the experimentally 
observed “dog house" and thinning of the coolant-channel 
wall could be captured with a single internal state variable 
model. Next, we will analyze the effects and advantages of 
reinforcing the cylindrical thrust chamber by placing the fibers 
along different orientations. Then we will try to find the 
optimum orientation of the fibers — the orientation that would 
yield the longest service life for the cylindrical thrust chambers. 

Conclusions 

A finite element viscoplastic stress-strain analysis for an 
experimental cyclindrical thrust chamber of NARloy-Z was 


made by using a viscoplastic model developed by Robinson. 
Computations were performed for two types of loading cycles. 
The short cycle of 3.5-sec duration pertained to the experi- 
mental situation, whereas the extended cycle of 485.1-sec 
duration corresponded to the operating cycle experienced by 
the space shuttle main engine. The analysis qualitatively 
replicated the “dog house" effect and the thinning of the 
coolant-channel wall for both the short and extended loading 
cycles. We noticed the “dog house" effect was caused by the 
biasing pressure differential across the coolant-channel wall. 
The location of extensive deformation in the component depen- 
ded on the duration of the loading cycle, thereby implying that 
failure of the component may also depend on the duration of 
the loading cycle. Futhermore, the results showed that 
structural analysis by a viscoplastic model that incorporates 
a single internal state variable, representing kinematic 
hardening, was capable of qualitatively predicting the experi- 
mentally observed behavior of the cylindrical thrust chamber— 
as was a viscoplastic model with two internal state variables 
(Freed's model). 
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